DailyExercise21_CacheLaPoudreData

Download Data

library(dataRetrieval)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(tidymodels)

── Attaching packages ────────────────────────────────────── tidymodels 1.2.0 ──
✔ broom        1.0.7     ✔ rsample      1.2.1
✔ dials        1.3.0     ✔ tune         1.2.1
✔ infer        1.0.7     ✔ workflows    1.1.4
✔ modeldata    1.4.0     ✔ workflowsets 1.1.0
✔ parsnip      1.2.1     ✔ yardstick    1.3.1
✔ recipes      1.1.0     
── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
✖ scales::discard() masks purrr::discard()
✖ dplyr::filter()   masks stats::filter()
✖ recipes::fixed()  masks stringr::fixed()
✖ dplyr::lag()      masks stats::lag()
✖ yardstick::spec() masks readr::spec()
✖ recipes::step()   masks stats::step()
• Search for functions across packages at https://www.tidymodels.org/find/

library(dplyr)
library(tsibble)

Registered S3 method overwritten by 'tsibble':
  method               from 
  as_tibble.grouped_df dplyr

Attaching package: 'tsibble'

The following object is masked from 'package:lubridate':

    interval

The following objects are masked from 'package:base':

    intersect, setdiff, union

# Example: Cache la Poudre River at Mouth (USGS site 06752260)
poudre_flow <- readNWISdv(siteNumber = "06752260",    
                          # Download data from USGS for site 06752260
                          parameterCd = "00060",      
                          # Parameter code 00060 = discharge in cfs)
                          startDate = "2013-01-01",   
                          # Set the start date
                          endDate = "2023-12-31") |>
  # Set the end date
  renameNWISColumns() |>                              
  # Rename columns to standard names (e.g., "Flow", "Date")
  mutate(Date = yearmonth(Date)) |>                   
  # Convert daily Date values into a year-month format (e.g., "2023 Jan")
  group_by(Date) |>                                   
  # Group the data by the new monthly Date
  summarise(Flow = mean(Flow))

GET:https://waterservices.usgs.gov/nwis/dv/?site=06752260&format=waterml%2C1.1&ParameterCd=00060&StatCd=00003&startDT=2013-01-01&endDT=2023-12-31

  # Calculate the average daily flow for each month

1. Convert to tsibble

# Convert to tsibble (monthly data, with Date as index)
poudre_ts <- poudre_flow |> 
  as_tsibble(index = Date)

2. Plotting the time series

library(ggplot2)
library(plotly)

Attaching package: 'plotly'

The following object is masked from 'package:ggplot2':

    last_plot

The following object is masked from 'package:stats':

    filter

The following object is masked from 'package:graphics':

    layout

# Basic ggplot of flow over time
p <- ggplot(poudre_ts, aes(x = Date, y = Flow)) +
  geom_line(color = "blue") +
  labs(title = "Monthly Average Streamflow - Cache la Poudre River",
       y = "Flow (cfs)", x = "Date") +
  theme_minimal()

# Animate with plotly
ggplotly(p)

3. Subseries

library(feasts)

Loading required package: fabletools

Attaching package: 'fabletools'

The following object is masked from 'package:yardstick':

    accuracy

The following object is masked from 'package:parsnip':

    null_model

The following objects are masked from 'package:infer':

    generate, hypothesize

# Subseries plot to show seasonal patterns
poudre_ts |> 
  gg_subseries(Flow)

4. Decompose

library(fabletools)

# STL decomposition with a seasonal window (say 13 months for annual seasonality)
decomp <- poudre_ts |> 
  model(STL(Flow ~ season(window = "periodic"))) |> 
  components()

# Plot the components
autoplot(decomp)